#!/usr/bin/octave -qp /home/trflp/web/lib/octave

# Durbin-Watson statistic
function s = dw(original, smoothed)
  n = length(original);
  d = original - smoothed;
  s = sum((d(2:end) - d(1:end-1)) .^ 2)/sum(d .^ 2);
  s = s*n/(n - 1);
endfunction

arg = argv;
sample = load("-ascii", arg{1});

				# calculate the optimal smoothing parameters by setting
				# the order to a fixed value and optimizing the window until
				# the Durbin-Watson statistic approaches the value of 2,
				# indicating zero correlation between adjacent data points
order = 2;
dwt = repmat(Inf, 1, 31);
# for window = 5:2:31
#   sm = sgolayfilt(sample, order, window);
#   dwt(window) = dw(sample, sm);
#   if dwt(window) < 1
#     break;
#   endif
# endfor
#[d, window] = min(abs(dwt - 2));
window = 9;
sm = sgolayfilt(sample, order, window); # redo smoothing with the optimal window

save -ascii - sm;
